สถิติแบบ Frequentist

สถิติแบบดั้งเดิมหรือที่เรียกกันว่า frequentist เป็นสถิติวิเคราะห์ที่ถูกใช้งานเป็นหลักในปัจจุบัน สถิติเกือบทุกตัวที่ผู้เรียนได้เรียนมาจนถึงปัจจุบันนี้ล้วนเป็นสถิติวิเคราะห์ภายใต้กระบวนทัศน์แบบ frequentist ทั้งสิ้น สถิติแบบ frequentist นี้มีข้อตกลงเบื้องต้นว่า (1) พารามิเตอร์เป็นค่าคงที่ ที่มีอยู่จริงในประชากรเป้าหมาย อย่างไรก็ตามด้วยข้อจำกัดที่ไม่สามารถเก็บรวบรวมข้อมูลจากทั้งประชากรได้ทำให้ผู้วิเคราะห์ไม่ทราบค่าพารามิเตอร์ดังกล่าว และจะประมาณโดยใช้ข้อมูลจากตัวอย่าง (2) ข้อมูลตัวอย่าง (sample data) เขียนแทนด้วย \(\bf{x}= \left \{\bf{x_1}, \bf{x_2},...,\bf{x_n} \right \}\) เป็นค่าที่สุ่มมาจากการประชากร ดังนั้นตัวอย่างที่สุ่มได้จึงเป็นค่าสุ่มที่และมีความเป็นไปได้ทั้งหมดเท่ากับ เมื่อ \(N\) คือขนาดประชากร และ \(n\) คือขนาดตัวอย่าง จากข้อตกลงเบื้องต้นข้อที่ (2) นี้จึงทำให้ (3) ค่าสถิติที่ประมาณได้จากข้อมูลตัวอย่างเป็นตัวแปรสุ่มที่มีการแจกแจงความน่าจะเป็น เรียกการแจกแจงความน่าจะเป็นของค่าสถิตินี้ว่า การแจกแจงความน่าจะเป็นของฟังก์ชันที่ได้จากตัวอย่างสุ่ม (sampling distribution)

ลักษณะการแจกแจงความน่าจะเป็นของค่าสถิติดังกล่าวมักประมาณได้ด้วย 2 วิธีการ วิธีการแรก ประมาณโดยใช้ทฤษฎีลิมิตลู่เข้าสู่ส่วนกลาง (central limit theorem) ซึ่งเป็นการคาดการณ์ลักษณะการแจกแจงของค่าสถิติในเชิงทฤษฎี เช่น การประมาณช่วงความเชื่อมั่นของค่าเฉลี่ยมีสูตรเป็น \(\overline{X}\pm t_{\alpha/2,n-1} \frac{S}{\sqrt{n}}\) จะเห็นว่าการประมาณช่วงความเชื่อมั่นดังกล่าวมีการใช้การแจกแจงที (student t distribution) ในการคำนวณส่วน margin of error การแจกแจงดังกล่าวพิสูจน์ในทางทฤษฎีโดยใช้ central limit theorem ได้ว่า เป็น sampling distribution ของ \(\overline{X}\) หากข้อมูลตัวอย่างที่ใช้ในการวิเคราะห์สุ่มมาจากประชากรที่มีการแจกแจงแบบปกติหรือใกล้เคียง แต่ไม่ทราบค่าความแปรปรวนของประชากร หรือถ้าหากประชากรไม่ได้มีการแจกแจงแบบปกติ แต่ตัวอย่างมีขนาดใหญ่เพียงพอ ก็ยังสามารถใช้ sampling distribution ดังกล่าวประมาณการแจกแจงของ \(\overline{X}\) ได้อยู่ หรือในการทดสอบสมมุติฐาน \(H_0: \mu=\mu_0\) สถิติทดสอบที่ใช้คือ \(t^*=\frac{\overline{X}-\mu_0}{S/ \sqrt{n}}\) จะมี sampling distributionในทางทฤษฎีเป็นการแจกแจงแบบที ที่มีองศาความเป็นอิสระเท่ากับ \(n-1\) เช่นเดียวกันหากข้อกำหนดเบื้องต้นของ central limit theorem เป็นจริงหรือใกล้เคียง การทราบการแจกแจงของค่าสถิติดังกล่าวทำให้ผู้วิเคราะห์สามารถประเมินความคลาดเคลื่อนในการตัดสินใจของการทดสอบสมมุติฐาน และสรุปผลได้ดังที่ได้เรียนไปแล้วในรายวิชาพื้นฐาน

ความถูกต้องของทฤษฎีนี้ขึ้นอยู่กับปัจจัยแวดล้อมที่เกี่ยวข้องว่ามีความสอดคล้องกับทฤษฎีมากน้อยเพียงใด ปัจจัยหนึ่งที่มีความสำคัญคือขนาดตัวอย่าง กล่าวคือหากขนาดตัวอย่างมีค่าน้อยเกินไป หรือการแจกแจงของข้อมูลมีความเบี่ยงเบนออกไปจากข้อตกลงเบื้องต้นของการวิเคราะห์ไปมาก ค่าสถิติที่คำนวณได้อาจไม่มีลักษณะการแจกแจงที่เป็นไปตามทฤษฎี และส่งผลให้การอนุมานเชิงสถิติมีความคลาดเคลื่อน อีกวิธีการหนึ่งที่สามารถใช้ในการประมาณ sampling distribution ของค่าสถิติได้คือใช้วิธีการสุ่มซ้ำ (resampling method) เช่น boostraping หรือ jackknife เพื่อประมาณแนวโน้มการแจกแจงของค่าสถิติ เป็นต้น

t.test(dat$Ach, mu=5.5)

## 
##  One Sample t-test
## 
## data:  dat$Ach
## t = 5.078, df = 149, p-value = 1.123e-06
## alternative hypothesis: true mean is not equal to 5.5
## 95 percent confidence interval:
##  5.709732 5.976934
## sample estimates:
## mean of x 
##  5.843333

คำถาม

1. output ข้างต้นแสดงผลการประมาณช่วงความเชื่อมั่น และทดสอบสมมุติฐานทางสถิติเกี่ยวกับผลสัมฤทธิ์ทางการเรียนของนักเรียน (คะแนนเต็ม 10) โดยจากการประมาณด้วยช่วงความเชื่อมั่น 95% ในข้างต้นพบว่ามีค่าเท่ากับ \([5.71, 5.98]\) ซึ่งหมายความว่าอะไร?

2. จากผลการทดสอบสมมุติฐาน \(H_0: \mu_{Ach}=5.5\) พบว่ามีค่า p-value < .000 หมายความว่าอะไร?

3. เราสามารถตัดสินใจยอมรับสมมุติฐาน \(H_0\) ได้หรือไม่?

จากข้อตกลงเบื้องต้นของสถิติแบบ frequentist นี้ทำให้การวิเคราะห์ข้อมูลเกิดข้อจำกัดต่าง ๆ ดังนี้

1. ด้านที่สำคัญคือด้านการแปลผลการวิเคราะห์ที่จะเห็นว่าความหมายของผลการวิเคราะห์ต่าง ๆ เป็นการกล่าวถึงพารามิเตอร์ที่สนใจแบบอ้อม ๆ ทั้งสิ้น ไม่มีการวิเคราะห์ใดที่สามารถอ้างอิงไปยังพารามิเตอร์ที่สนใจได้โดยตรง ทั้งนี้เป็นเพราะการอนุมานเชิงสถิติแบบ frequentist พึ่งพาเครื่องมือหลักคือ samping distribution ของค่าสถิติ ในขณะที่พารามิเตอร์ถูกสมมุติให้เป็นค่าคงที่ตั้งแต่ต้น probability statement ต่าง ๆ จึงเป็นของค่าสถิติทั้งสิ้น ดังความความเชื่อมั่นที่รับประกันในช่วงความเชือมั่น และค่าความคลาดเคลื่อนในการทดสอบสมมุติฐานจึงเป็นค่าความน่าจะเป็นที่ใช้รับประกันความเป็นไปได้ของค่าสถิติทั้งสิ้นไม่ใช่ของค่าพารามิเตอร์ที่สนใจจริง ๆ ในการวิเคราะห์

2. ด้านที่สองคือด้านการประมาณค่าพารามิเตอร์ภายใต้โมเดลซับซ้อน เครื่องมือสำคัญสำหรับประมาณค่าพารามิเตอร์ในโมเดลต่าง ๆ ของสถิติแบบ frequentist คือตัวประมาณ เช่น ตัวประมาณกำลังสองน้อยสุด ตัวประมาณภาวะความควรจะเป็นสูงสุด ซึ่งล้วนเป็นวิธีการที่จะหาค่าประมาณที่ดีที่สุดโดยอิงจากจุดต่ำสุด หรือจุดสูงสุดของฟังก์ชันวัตถุประสงค์ ดังตัวอย่างในรูปด้านล่าง ซึ่งในกรณีที่โมเดลมีความซับซ้อนมาก ๆ ฟังก์ชันวัตถุประสงค์ดังกล่าวจะมีความซับซ้อนโดยมีจำนวนมิติที่สูงขึ้น และอาจมีจุดสูงต่ำ และ/หรือ จุดต่ำสุดสัมพัทธ์จำนวนมาก ซึ่งเป็นอุปสรรคให้การประมาณค่าด้วยวิธีการในข้างต้นอาจได้ค่าประมาณที่คลาดเคลื่อนไปจากค่าที่เหมาะสมที่สุด หรืออัลกอริทึมไม่สามารถหาค่าประมาณที่เหมาะสมได้ภายใต้จำนวนรอบของการประมาณที่กำหนด

3. สืบเนื่องจากที่การอนุมานเชิงสถิติแบบ frequentist นั้นมักอิงกับทฤษฎีความน่าจะเป็นที่อิงข้อกำหนดเบื้องต้นที่ขนาดตัวอย่างต้องมีขนาดใหญ่เพียงพอ ซึ่งทำให้สถิติแบบ frequentist มักทำงานได้ไม่ดีในกรณีที่ตัวอย่างมีขนาดเล็ก

สถิติแบบเบส์ (Bayesian Statistics)

สถิติแบบเบส์เป็นกระบวนทัศน์ในการวิเคราะห์ข้อมูลอีกลักษณะหนึ่งที่มีความแตกต่างจากไปกระบวนทัศน์แบบดั้งเดิมที่นักวิเคราะห์ข้อมูลทั่วไปคุ้นเคยกัน ซึ่งทำให้สถิติแบบเบส์มีจุดเด่นและสามารถแก้ไขข้อจำกัดที่พบในสถิติแบบ frequentist ได้ ความแตกต่างดังกล่าวเริ่มตั้งแต่ข้อตกลงเบื้องต้นของสถิติแบบเบส์ที่แตกต่างจากสถิติแบบดั้งเดิมอย่างมาก กล่าวคือ

(1) พารามิเตอร์เป็นตัวแปรสุ่ม สถิติแบบเบส์มองว่าพารามิเตอร์เป็นค่าที่ผู้วิเคราะห์สนใจจะศึกษาและต้องการทราบอย่างแท้จริง แต่อย่างไรก็ตามเนื่องจากผู้วิเคราะห์ไม่ทราบข้อมูลประชากร การคาดการณ์หรือความเชื่อ (belief) เกี่ยวกับพารามิเตอร์ดังกล่าวจึงมีความไม่แน่นอน สภาพดังกล่าวจึงสมเหตุสมผลที่จะกำหนด probability statement ให้กับพารามิเตอร์ในโมเดล

จากข้อกำหนดนี้ผลการประมาณค่าพารามิเตอร์ที่ได้จากวิธีการแบบเบส์จึงจะไม่ได้ค่าประมาณแบบค่าเดียว แต่จะได้เป็นการประมาณการแจกแจงของพารามิเตอร์แทน เช่น การประมาณค่าเฉลี่ยผลสัมฤทธิ์ของนักเรียน ผลการวิเคราะห์แบบเบส์จะไม่ได้ให้เป็นค่าประมาณเพียงค่าเดียวคือ 5.84 คะแนน ดัง output ข้างต้น แต่จะให้เป็นการแจกแจงของพารามิเตอร์ค่าเฉลี่ยผลสัมฤทธิ์ทางการเรียน ซึ่งเรียกว่าการแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) และสารสนเทศเกี่ยวกับค่าเฉลี่ยของผลสัมฤทธิ์ดังกล่าวจะถูกสกัดออกมาจากการแจกแจงนี้

(2) ข้อมูลตัวอย่างเป็นค่าคงที่ ถึงแม้ว่าในความเป็นจริงตัวอย่างที่สุ่มมาจากประชากรจะมีความเป็นไปได้ที่หลากหลาย และมีความไม่แน่นอนในการได้มาซึ่งตัวอย่างแต่ละชุด อย่างไรก็ตามเมื่อดำเนินการสุ่มตัวอย่างมาทำการวิเคราะห์แล้ว ความเป็นไปได้จำนวนมากก่อนการสุ่มตัวอย่างจะเหลือเพียงความเป็นไปได้เดียว สถิติแบบเบส์จึงถือว่าข้อมูลตัวอย่างเป็นค่าคงที่ ไม่มีการกำหนด probability statetment ให้กับข้อมูลตัวอย่างนี้

สถิติแบบเบส์สามารถนำมาใช้เพื่อแก้ไขข้อจำกัดที่เกิดขึ้นจากสถิติแบบ frequentist ได้ดังนี้

1. จากข้อกำหนดเบื้องต้นข้างต้น การอนุมานเชิงสถิติด้วยสถิติแบบเบส์จะดำเนินการผ่านการแจกแจงความน่าจะเป็นของพารามิเตอร์โดยตรงที่เรียกว่า การแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) การแจกแจงดังกล่าวเป็นเครื่องมือสำคัญในการอนุมานเชิงสถิติแบบเบส์ แทน sampling distribution ในการอนุมานเชิงสถิติแบบดั้งเดิม ซึ่งจุดนี้ทำให้สถิติแบบเบส์สร้างความแตกต่างอย่างมากเมื่อเปรียบเทียบกันสถิติแบบดั้งเดิม โดยจะเห็นว่าการประมาณค่าพารามิเตอร์แบบเบส์ไม่ได้อิงกับค่าสถิติเพียงค่าเดียวเหมือนกับสถิติแบบดั้งเดิม แต่อิงกับการแจกแจงความน่าจะเป็นภายหลัง

2. การที่สถิติแบบเบส์ใช้การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์เป็นเครื่องมือหลักในการอนุมานเชิงสถิติ จึงทำให้การวิเคราะห์ทำได้อย่างตรงไปตรงมาโดยไม่ต้องอาศัยทฤษฎีหรือวิธีการที่ซับซ้อนเหมือนกับสถิติแบบดั้งเดิม กล่าวคือผู้วิเคราะห์ใช้เพียงสถิติพื้นฐาน เช่น ค่าเฉลี่ย มัธยฐาน ฐานนิยม ส่วนเบี่ยงเบนมาตรฐาน หรือเปอร์เซ็นไทล์ เป็นเครื่องมือในการสรุปสารสนเทศเกี่ยวกับพารามิเตอร์ที่สนใจออกมาจากการแจกแจงความน่าจะเป็นภายหลังที่ประมาณได้ ก็เพียงพอแล้ว นอกจากนี้การแปลผลเพื่ออนุมานเกี่ยวกับพารามิเตอร์ยังทำได้ง่ายและตรงไปตรงมา และได้ข้อสรุปโดยตรงไปยังปริภูมิของพารามิเตอร์ที่เป็นเป้าหมาย แตกต่างจากข้อสรุปที่ได้จากสถิติแบบดั้งเดิมที่เป็นข้อสรุปเกี่ยวกับการแจกแจงของค่าสถิติซึ่งอยู่บนปริภูมิตัวอย่างเท่านั้น

3. การประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ในโมเดล จะใช้สารสนเทศสองส่วนด้วยกัน ส่วนแรกเป็นสมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ที่สนใจ ซึ่งกำหนดอยู่ในรูปการแจกแจงความน่าจะเป็นที่เรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และส่วนที่สองคือสารสนเทศที่ได้จากข้อมูลตัวอย่าง ซึ่งแตกต่างจากสถิติแบบดั้งเดิมที่จะใช้สารสนเทศจากข้อมูลตัวอย่างเท่านั้น

ลองพิจารณาตัวอย่างต่อไปนี้

สมมุติว่าต้องการคาดการณ์คะแนนสอบของนายบูม ที่มีความเป็นไปได้ 4 แบบดังรูป

ก่อนการเก็บรวบรวมข้อมูลจริงในชั้นเรียนของนายบูม สมมุติว่าผู้วิเคราะห์มีข้อมูลในอดีตที่เกี่ยวข้องกับคะแนนสอบของนายบูมดังนี้

1. ไม่มีข้อมูลใด ๆ ที่เป็นประโยชน์ (ให้น้ำหนักกับความเป็นไปได้ A1, A2, A3 และ A4 อย่างไร?)

2. หากทราบเพิ่มเติมว่าการสอบที่ผ่าน ๆ มาโดยเฉลี่ยนายบูมสอบได้คะแนนคิดเป็นอันดับที่ 4 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 จะเหมือนเดิมมั้ย?)

3. หากทราบเพิ่มเติมอีกว่า การสอบที่ผ่าน ๆ มา แก้วซึ่งเป็นเพื่อนในห้องเรียนของนายบูม สอบได้คะแนน 38 คะแนน คิดเป็นคะแนนที่น้อยที่สุดในชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

4. หากทราบอีกว่า ในการสอบที่ผ่าน ๆ มา นิดสอบได้คะแนนโดยเฉลี่ยคิดเป็น 89 คะแนน ซึ่งเป็นอันดับที่ 3 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

การแจกแจงก่อนหน้าของคะแนนสอบนายบูม

ผู้วิเคราะห์สามารถปรับสารสนเทศในการแจกแจงความน่าจะเป็นก่อนหน้าในข้างต้นได้ หากมีการเก็บรวบรวมข้อมูลตัวอย่างเพิ่มเติม เช่น

สมมุติว่าการแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมเป็นแบบ uniform และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

อีกตัวอย่างหนึ่ง สมมุติว่า การแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมมีน้ำหนักสูงที่ A3 และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

กำหนดตัวแปรในการวิเคราะห์

สมมุติว่าผู้วิจัยต้องการวิเคราะห์ความลำเอียงในการออกหน้าหัวหรือก้อยของเหรียญที่ผลิตขึ้นมารุ่นหนึ่ง การวิเคราะห์นี้ตัวแปรที่สนใจจึงเป็นผลลัพธ์ที่สังเกตได้จากการโยนเหรียญในแต่ละครั้ง กำหนดให้ \(y_i\) แทนค่าสังเกตที่ได้จากการโยนเหรียญในครั้งที่ \(i\) จะได้ว่า

\(y_i=\left\{\begin{matrix} 1 & ; H\\ 0 & ; T \end{matrix}\right.\)

กำหนดโมเดลของค่าสังเกต

เมื่อกำหนดตัวแปรที่สนใจแล้วจะพบว่าตัวแปรที่สนใจมีเพียงตัวแปรตาม 1 ตัว โดยมีค่าสังเกตที่เป็นไปได้เพียง 2 ค่าได้แก่ H หรือ T ซึ่งแทนด้วย 1 กับ 0 ตามลำดับ โมเดลของค่าสังเกต \(y_i\) นี้จึงสามารถกำหนดได้ด้วยการแจกแจงความน่าจะเป็นแบบ Bernoulli กล่าวคือ ความน่าจะเป็นที่จะออกหน้าหัวมีค่าเท่ากับ \(p(y_i=1|\theta)=\theta\) และความน่าจะเป็นที่ออกหน้าก้อยมีค่าเท่ากับ \(p(y_i=0|\theta)=1-\theta\) ซึ่งเขียนเป็นสมการรวมได้ดังนี้

\(p(y_i|\theta)=\theta^{y_i}(1-\theta)^{1-y_i}\) เมื่อ \(i = 1,2,3,...,n\)

จากโมเดลในข้างต้นจะเห็นว่าการเกิดค่าสังเกต \(y_i\) ขึ้นอยู่กับพารามิเตอร์ \(\theta\) ที่มีค่าอยู่บนช่วง \([0,1]\) หากพารามิเตอร์ดังกล่าวมีค่าเท่ากับ \(\theta=0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวและก้อยมีค่าเท่ากัน แต่ถ้าหากพารามิเตอร์ \(\theta>0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวมีมากกว่า จากความหมายดังกล่าวทำให้สามารถมองได้ว่าพารามิเตอร์ \(\theta\) เป็นพารามิเตอร์ความลำเอียงของเหรียญ หากประมาณพารามิเตอร์ดังกล่าวได้จะสามารถสรุปความลำเอียงของเหรียญได้

เนื่องจากเราไม่ได้เก็บรวบรวมข้อมูลเพียงค่าสังเกตเดียว แต่เราทำการทดลองซ้ำ ๆ จำนวน \(n\) ครั้ง เนื่องจากการทดลองแต่ละครั้งเป็นอิสระซึ่งกันและกัน ทำให้โมเดลของชุดค่าสังเกตหรือเรียกอย่างเป็นทางการว่า ฟังก์ชันภาวะความควรจะเป็นดังนี้

\(p(\bf{y}|\theta)=p(y_1|\theta)p(y_2|\theta)...p(y_n|\theta)\)

ซึ่งมีค่าเท่ากับ

\(p(\bf{y}|\theta)=\theta^{\sum_{i=1}^ny_i}(1-\theta)^{n-\sum_{i=1}^ny_i}\)

ในฟังก์ชันภาวะความควรจะเป็นนี้จะให้ค่าความน่าจะเป็นหรือค่าความหนาแน่นของข้อมูลตัวอย่าง \(\bf{y}\) บนแต่ละค่าที่เป็นไปได้ของพารามิเตอร์ \(\theta\) เนื่องจากข้อมูลตัวอย่างเป็นค่าคงที่ ดังนั้นค่าความหนาแน่นนี้จึงเปลี่ยนแปลงไปตามค่าของพารามิเตอร์ \(\theta\) โดยหากค่าความหนาแน่นนี้มีค่าสูง ณ ค่า \(\theta\) ค่าใด นั่นหมายถึงว่าพารามิเตอร์ค่านั้นทำให้โมเดลของค่าสังเกตกับค่าสังเกตจริงมีความสอดคล้องกันสูง ข้อสังเกตนี้ถูกนำไปใช้เพื่อหาค่าประมาณด้วยวิธีภาวะความควรจะเป็นสูงสุด (maximum likelihood: ML) แต่ในวิธีการแบบเบส์จะนำสารสนเทศส่วนนี้ไปประมวลร่วมกับความเป็นไปได้ของพารามิเตอร์ในการแจกแจงความน่าจะเป็นเบื้องต้นก่อน

กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์

การกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) สามารถทำได้หลายลักษณะ ทั้งแบบไม่ต่อเนื่อง และแบบต่อเนื่อง ขึ้นอยู่กับขอบเขตของการวิเคราะห์ ข้อมูลในอดีตที่ใช้สนับสนุน และธรรมชาติของพารามิเตอร์นั้น ในตัวอย่างข้างต้นพารามิเตอร์ \(\theta\) ต้องมีค่าอยู่บนช่วง \([0,1]\)

เพื่อให้ตัวอย่างนี้ทำความเข้าใจได้ง่าย จึงเลือกกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) เป็นการแจกแจงแบบไม่ต่อเนื่อง โดยมีค่าที่เป็นไปได้เท่ากับ 0.0, 0.1, 0.2, 0.3, …., 1.0 และมีการแจกแจงของความน่าจะเป็นก่อนหน้าดังรูป

# prior distribution
theta<-seq(0,1,0.1)
p.theta<-c(0.01,0.04,0.075,0.1,0.15,0.25,0.15,0.1,0.075,0.04,0.01)
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab=expression(p(theta)), 
    type="h", lwd=2, col="orange",
    main="Prior Distribution")
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

การประมาณการแจกแจงความน่าจะเป็นภายหลัง

จากการทดลองโยนเหรียญจำนวน 10 ครั้งพบว่า ออกหน้าหัวจำนวน 6 ครั้ง ดังนั้นจะได้ว่า

\(p(\bf{y}|\theta)=\theta^{6}(1-\theta)^4\)

จากขอบเขตของ \(theta\) ที่กำหนดในการแจกแจงความน่าจะเป็นเบื้องต้นจะได้ว่าฟังก์ชันภาวะความควรจะเป็น มีลักษณะดังรูป

# likelihood function
L<-theta^6*(1-theta)^4
plot(theta,L, 
     xlab=expression(theta), 
    ylab="p(y|theta)",
    type="h", lwd=2, col="orange",
    main="Likelihood Function")
points(theta,L, type="p",pch=16, cex=2, col="orange")

จากทฤษฎีบทของเบส์เราสามารถรวมสารสนเทศจากการแจกแจงความน่าจะเป็นเบื้องต้น กับฟังก์ชันภาวะความควรจะเป็นได้ดังนี้

# marginal y
p.y<-sum(p.theta*L)

par(mfrow=c(1,2))
# posterior of theta given y
post.theta<-L*p.theta/p.y
plot(theta,post.theta, 
     xlab=expression(theta), 
    ylab="p(theta|y)",
    type="h", lwd=2, col="blue",
    ylim=c(0,0.5),
    main="Posterior Distribution")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")

# prior and posterior in the same plot
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab="Probability", 
    type="l", lwd=2, col="orange",
    ylim=c(0,0.5))
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

points(theta,post.theta, 
    type="l", lwd=2, col="blue")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
legend(-0.05,0.48, legend=c("Prior","Posterior"), col=c("orange","blue"), lty=1, bty="n", lwd=2)

การอนุมานเชิงสถิติจากการแจกแจงความน่าจะเป็นภายหลัง

เมื่อได้การแจกแจงความน่าจะเป็นภายหลังที่ต้องการแล้ว ผู้วิเคราะห์สามารถวิเคราะห์การแจกแจงดังกล่าวเพื่ออนุมานเกี่ยวกับพารามิเตอร์ความลำเอียงที่ต้องการศึกษาได้ โดยสามารถทำได้หลายลักษณะ ทั้งการประมาณค่าแบบจุด การประมาณค่าแบบช่วง และการทดสอบสมมุติฐาน รายละเอียดดังนี้

(1) การประมาณค่าแบบจุด

การประมาณค่าแบบจุดสามารถทำได้โดยใช้สถิติพื้นฐานสรุปแนวโน้มสู่ส่วนกลางของค่าพารามิเตอร์ จากการแจกแจงความน่าจะเป็นภายหลัง เช่น ค่าเฉลี่ย มัธยฐาน และฐานนิยม นอกจากนี้ยังสามารถคำนวณค่าส่วนเบี่ยงเบนมาตรฐานเพื่อประเมินระดับความน่าเชื่อถือของค่าประมาณแบบจุดดังกล่าวได้อีกด้วย จากการแจกแจงความน่าจะเป็นภายหลังข้างต้น จะได้ว่า

# mean = expected value of theta
mean.theta<-sum(post.theta*theta)
sd.theta<-sqrt(sum(post.theta*(theta-mean.theta)^2))
mean.theta

## [1] 0.5540441

sd.theta

## [1] 0.1145753

(2) การประมาณค่าแบบช่วง

ช่วงการประมาณที่ได้จากวิธีแบบเบส์จะเรียกว่า ช่วงความน่าเชื่อถือ (credible interval) การประมาณค่าพารามิเตอร์แบบช่วงจากการแจกแจงความน่าจะเป็นภายหลังสามารถทำได้หลายวิธีการ วิธีการง่าย ๆ คือการใช้ค่าเฉลี่ยบวกลบด้วยส่วนเบี่ยงเบนมาตรฐานของพารามิเตอร์ อีกวิธีการหนึ่งที่มีประสิทธิภาพมากกว่าคือการหาช่วงแบบ highest density interval (HDI) กล่าวคือเป็นช่วงการประมาณที่มีความหนาแน่นมากที่สุดที่ทำให้ค่าความน่าจะเป็นของพารามิเตอร์ภายในช่วงดังกล่าวเท่ากับค่าที่กำหนดเช่น ช่วง 95% HDI ของพารามิเตอร์ \(\theta\) ในข้างต้นจะมีค่าประมาณ \([0.4,0.8]\)

ช่วงดังกล่าวมีความแตกต่างจาก 95% confidence interval อย่างไร??

(3) การทดสอบสมมุติฐาน

การทดสอบสมมุติฐานด้วยวิธีการแบบเบส์มีความยืดหยุ่นสูงมาก และสามารถทำได้หลายวิธีการ วิธีการแรกเรียกว่า maximum a posteriori (MAP) test มีรายละเอียดดังนี้

กำหนดให้ \(H_0\) และ \(H_1\) เป็นสมมุติฐานที่ต้องการเปรียบเทียบ ในกรณีนี้จะเห็นว่ามีความไม่แน่นอนว่าสมมุติฐานใดเป็นสมมุติฐานที่ถูกต้องกันแน่ ดังนั้นจึงสามารถกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานทั้งสองได้ ดังนี้

\(p(H_0) = p_0\) และ \(p(H_1)=p_1\) โดยที่ \(p_0+p_1=1\)

และกำหนดให้ \(p(\bf{y}|H_0)\) กับ \(p(\bf{y}|H_1)\) คือฟังก์ชันภาวะความควรจะเป็นบนสมมุติฐานทั้งสอง ซึ่งจากทฤษฎีบทของเบส์จะได้ว่า

\(p(H_0|\bf{y})=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y})}\)

\(p(H_1|\bf{y})=\frac{p(\bf{y}|H_1)p(H_1)}{p(\bf{y})}\)

การตัดสินใจว่าควรเลือกเชื่อสมมุติฐานใด สามารถทำได้โดยดูจากอัตราส่วนที่เรียกว่า posterior odd ซึ่งมีค่าเท่ากับ

\(\frac{p(H_0|\bf{y})}{p(H_1|\bf{y})}=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y}|H_1)p(H_1)}\)

โดยหากอัตราส่วนข้างต้นมีค่ามากกว่า 1 แสดงว่า สมมุติฐาน \(H_0\) มีแนวโน้มน่าเชื่อถือมากกว่า ในทางกลับกันหากอัตราส่วนดังกล่าวมีค่าน้อยกว่า 1 แสดงว่าสมมุติฐาน \(H_1\) มีแนวโน้มน่าเชื่อถือมากกว่า

จากตัวอย่างโยนเหรียญ หากต้องการทดสอบว่าเหรียญมีความเที่ยงตรงหรือไม่ อาจกำหนดสมมุติฐานเป็นดังนี้

\(H_0: \theta = 0.5\) vs \(H_1: \theta>0.5\)

prior.H<-c(0.5,0.5) # prior for H0 and H1 respectively.
# calculate likelihood
L.H0<-0.5^6*(1-0.5)^4
theta.H1<-theta[theta>0.5]
L.H1<-sum(theta.H1^6*(1-theta.H1)^4)

# calculate posterior odd
(prior.H[1]*L.H0)/(prior.H[2]*L.H1)

## [1] 0.3727444

อีกวิธีการหนึ่งที่สามารถทำได้เรียกว่า minimum cost hypothesis test รายละเอียดมีดังนี้

กำหนดให้

• \(C_{10}\) คือ cost ของการเลือก \(H1\) โดยที่ \(H_0\) ถูกต้อง (cost ของการเกิด type I error)
• \(C_{01}\) คือ cost ของการเลือก \(H_0\) โดยที่ \(H_1\) ถูกต้อง (cost ของการเกิด type II error)

เกณฑ์การพิจารณาจะใช้ค่าความเสี่ยงภายหลัง (posterior risk) ดังนี้

\(\frac{p(H_0|\bf{y})C_{10}}{p(H_1|\bf{y})C_{01}}\)

จากตัวอย่างโยนเหรียญสมมุติว่า \(C_{10}=5C_{01}\)

c10<-5
c01<-1

# posterior risk of accepting H0
(prior.H[2]*L.H1)*c01

## [1] 0.001309962

# posterior risk of accepting H1
(prior.H[1]*L.H0)*c10

## [1] 0.002441406

# risk ratio
(prior.H[1]*L.H0)*c10/(prior.H[2]*L.H1)*c01

## [1] 1.863722

อีกวิธีการที่ง่ายกว่าสองวิธีการแรกมากคือ วิธี ROPE ซึ่งเป็นการเปรียบเทียบกันระหว่างช่วงความน่าเชื่อถือแบบ HDI กับช่วง ROPE ที่เรียกว่า region of practical equivalence ซึ่งเป็นช่วงที่ยอมรับได้ในทางปฏิบัติ (โดยผู้วิเคราะห์) ว่าหากพารามิเตอร์ยังมีค่าอยู่ในช่วง ROPE ดังกล่าวอยู่ แสดงว่ายังไม่แตกต่างจาก \(H_0\) โดยหากช่วง HDI ครอบคลุมช่วง ROPE เราจะสามารถยอมรับ \(H_0\) ได้ แต่หากช่วง HDI ไม่มีส่วนใดที่ครอบคลุมช่วง ROPE เลย แสดงว่าไม่ช่วงเชื่อถือ \(H_0\) อย่างยิ่ง

การคัดเลือกโมเดล (model selection)

การคัดเลือกโมเดลเป็นวิธีการหนึ่งที่สามารถใช้ในการทดสอบมสมมุติฐานของผู้วิเคราะห์ รวมทั้งใช้เปรียบเทียบและคัดเลือกโมเดลที่เหมาะสมที่สุดสำหรับอธิบายข้อมูล สถิติตัวหนึ่งที่มีความสำคัญคือ bayes factor วัตถุประสงค์ของ bayes factor ถูกใช้เพื่อวัดว่าข้อมูลตัวอย่างที่เก็บรวบรวมมมานี้สนับสนุนโมเดลหรือสมมุติฐานตัวใดมากกว่ากัน bayes factor มีค่าเท่ากับอัตราส่วนระหว่างค่าภาวะความควรจะเป็นของสมมุติฐานสองตัวบนชุดข้อมูลตัวอย่างเดียวกัน สมมุติให้ \(H_1\) คือโมเดลตามสมมุติฐานที่ 1 และ \(H_2\) คือโมเดลตามสมมุติฐานที่ 2 และ \(\bf{x}\) คือชุดข้อมูลตัวอย่าง จะได้ว่าสถิติ bayes factor (BF) มีค่าเท่ากับ

\(BF = \frac{p(\bf{x}|H_1)}{p(\bf{x}|H_2)}=\frac{p(H_1|\bf{x})}{p(H_2|\bf{x})}\times \frac{p(H_2)}{p(H_1)}\)

จากนิยามข้างต้นจะเห็นว่า bayes factor ตามนิยามดังกล่าวสามารถเขียนใหม่ในรูปผลคูณระหว่างอัตราส่วนสองตัว ได้แก่ อัตราส่วนระหว่างการแจกแจงความน่าจะเป็นภายหลังกับอัตราส่วนระหว่างการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานที่ 1 ต่อสมมุติฐานที่ 2

จากนิยามนี้สามารถแปลผลได้ว่าหาก \(BF\) มีค่ามากกว่า 1 นั่นหมายถึงข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_1\) มากกว่า ในทางกลับกันหาก \(BF\) มีค่าน้อยกว่า 1 นั่นหมายความว่าข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_2\) มากกว่า Harold Jeffery (1998) ได้เสนอสเกลการแปลความหมายสถิติ bayes factor ไว้ดังนี้

การคำนวณและการใช้งาน bayes factor จะกล่าวถึงในสัปดาห์ที่สองของหัวข้อนี้

นอกจากสถิติ bayes factor แล้วยังมีสถิติอีกหลายตัวที่สามารถใช้เปรียบเทียบและคัดเลือกโมเดลแบบเบส์ได้ เช่น deviance information criterion (DIC), widely applicable information criterion (WAIC) และ leave-one-out cross-validation (LOO) สถิติทั้งหมดที่กล่าวมานี้ใช้สะท้อนความแตกต่างกันระหว่างข้อมูลตัวอย่าง \(y\) กับค่าทำนายที่ generate จากโมเดลการวิเคราะห์ \(\hat{y}\) ค่าที่ได้จากสถิติดังกล่าวจึงใช้สะท้อนความคลาดเคลื่อนของแต่ละโมเดลเมื่อเปรียบเทียบกับข้อมูลจริง โมเดลที่มีค่าสถิติดังกล่าวต่ำกว่าจะเป็นโมเดลที่มีความเหมาะสมมากกว่า

การแจกแจงความน่าจะเป็นก่อนหน้า

การแจกแจงความน่าจะเป็นก่อนหน้าใช้สะท้อนความเชื่อเบื้องต้นในพารามิเตอร์ต่าง ๆ ของโมเดล ความเชื่อดังกล่าวมาได้จากทั้งผลการศึกษาในอดีต หลักเหตุผล และประสบการณ์ของนักวิจัย นอกจากนี้ยังควรต้องพิจารณาธรรมชาติหรือค่าที่เป็นไปได้ของพารามิเตอร์ประกอบด้วย

โดยทั่วไปการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าสามารถทำได้ 2 ลักษณะ ลักษณะแรกเรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ (non-informative prior distribution) และการแจกแจงความน่าจะเป็นก่อนหน้าแบบให้สารสนเทศ (informative prior distribution)

การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์อิทธิพลคงที่สามารถเลือกใช้ได้หลายการแจกแจง เช่น การแจกแจงแบบ uniform ที่เป็น noninformative prior การแจกแจงแบบปกติ ที่เป็นไปได้ทั้ง noninformative และ informative prior ขึ้นอยู่กับการกำหนดพารามิเตอร์ค่าเฉลี่ย (mean) และค่าความถูกต้อง (precision) โดยที่ค่าความถูกต้องดังกล่าวเป็นส่วนกลับของพารามิเตอร์ความแปรปรวน (variance) ดังนี้ precision = 1/variance

จากตัวอย่างข้างต้นกำหนดให้

\(\beta_0 \sim N(0,10^{-4})\)

\(\beta_1 \sim N(0,10^{-4})\)

การกำหนดข้างต้นคือการแจกแจงความน่าจะเป็นแบบปกติที่มีค่าเฉลี่ยเท่ากับ 0 และความแปรปรวนเท่ากับ 10000 ซึ่งเรียกได้ว่าเป็นการแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ รูปด้่านล่างแสดงลักษณะการแจกแจงแบบปกติข้างต้น จากรูปจะเห็นว่าการแจกแจงดังกล่าวมีขอบเขตที่กว้างพารามิเตอร์ส่วนใหญ่จะมีค่าตกอยู่ในช่วง -200 ถึง 200

การแจกแจงความน่าจะเป็นก่อนหน้าของส่วนเบี่ยงเบนมาตรฐานหรือความแปรปรวน มีความจำเป็นต้องเลือกการแจกแจงที่ domain ของการแจกแจงมีค่าไม่ติดลบ โดยปกติมักกำหนดให้ความแปรปรวนดังกล่าวมีการแจกแจงแบบแกมมาผกผัน (inverse gamma distribution)

จากตัวอย่างข้างต้น กำหนดให้

\(sigma^2 \sim IG(a,b)\)

โดยที่ \(IG\) คือการแจกแจงแบบ inverse-gamma

Note: ในเชิงทฤษฎีหากกำหนดให้ \(sigma^2 \sim IG(a,b)\) จะได้ว่า

\(p(\sigma^2|a,b) \propto (\sigma^2)^{-a-1}exp(\frac{-b}{\sigma^2})\)

หากกำหนดให้ \(a \rightarrow 0\) และ \(b \rightarrow 0\) แล้วฟังก์ชันความหนาแน่นในข้างต้นจะเขียนใหม่ได้เป็น

\(p(\sigma^2|a,b) \propto \frac{1}{\sigma^2}\) เรียกว่า Jeffery prior ซึ่งเป็น noninformative prior ที่ใช้ได้ตัวหนึ่งของพารามิเตอร์ความแปรปรวน

ดังนั้นในกรณีนี้จึงกำหนดให้

\(sigma^2 \sim IG(0.0001,0.0001)\)

รูปด้านล่างแสดงรายการของการแจกแจงความน่าจะเป็นที่สามารถกำหนดได้ใน BUGS และ JAGS

จากการระบุโมเดลข้างต้นทำให้สามารถเขียน syntax เพื่อระบุโมเดลใน JAGS ได้ดังนี้

model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~gamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}

การประมวลผลด้วย JAGS

การประมวลผลด้วย JAGS มี 4 ขั้นตอนหลักได้แก่

(1) เตรียมข้อมูล

(2) ระบุโมเดล

(3) ส่งข้อมูลและโมเดลไปประมวลผลบน JAGS

(4) ตรวจสอบการลู่เข้าของ MCMC

(5) วิเคราะห์และแปลผล

modelString<-"model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1*x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~dgamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}" #Wrap model syntax into String object

writeLines(modelString, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/simReg.txt") # write modelString to model.txt

ATT<-rnorm(100,30,10)
ACH<-30+0.55*ATT+rnorm(100,0,5)
dat.reg<-data.frame(ATT,ACH)

head(dat.reg)

##         ATT      ACH
## 1 32.133899 52.34851
## 2  8.208871 37.08404
## 3 37.542008 54.13731
## 4 13.544627 38.93896
## 5 35.487326 44.67005
## 6 42.563916 50.66076

par(mar=c(4,4,1,1))
plot(ATT,ACH, pch=16, xlab="ATT", ylab="ACH")

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(rjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model

update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 10000,
                      thin=1)
head(samples)

## [[1]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 30.78025 0.5202950 21.90599
## [2,] 30.60087 0.4840719 24.14094
## [3,] 31.66543 0.4708754 37.18908
## [4,] 31.32383 0.4505433 21.38748
## [5,] 31.89650 0.4600205 24.17811
## [6,] 31.78018 0.4470285 30.98406
## [7,] 32.44455 0.4417329 24.82733
## 
## [[2]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 29.98931 0.4762332 30.62127
## [2,] 29.89494 0.5395243 21.85398
## [3,] 29.32648 0.5483189 30.41757
## [4,] 29.36315 0.5429597 31.03104
## [5,] 29.14981 0.5379856 22.40325
## [6,] 30.14993 0.5327550 26.61984
## [7,] 29.80438 0.5197910 23.47122
## 
## [[3]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 31.81481 0.4581798 24.36047
## [2,] 31.77167 0.4908677 20.95915
## [3,] 31.31646 0.4657395 25.43489
## [4,] 31.52530 0.4829907 30.82938
## [5,] 31.05918 0.4862660 22.68679
## [6,] 30.97493 0.5190658 24.36354
## [7,] 30.30192 0.5135753 24.93864
## 
## attr(,"class")
## [1] "mcmc.list"

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(runjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

runjag.output<-run.jags(method="parallel",
                model="simReg.txt",
                monitor=c("b0","b1","sigma2"),
                data=dataList,
                inits=initsList,
                n.chains=3,
                burnin=2000,
                sample=10000, #n.iter
                thin=1,
                summarise=FALSE,
                plots=FALSE)
  
samples<-as.mcmc.list(runjag.output)
head(samples)
plot(samples)

library(MCMCvis)
MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

MCMCsummary(samples)

##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     30.0629002 1.61239595 26.9107602 30.0698484 33.2015106    1  1633
## b1      0.5213053 0.05160855  0.4217729  0.5211085  0.6222152    1  1607
## sigma2 26.9798622 3.93882192 20.3574546 26.6128567 35.6623178    1 27719

summary(samples)

## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## b0     30.0629 1.61240 0.009309       0.039919
## b1      0.5213 0.05161 0.000298       0.001288
## sigma2 26.9799 3.93882 0.022741       0.023665
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%    75%   97.5%
## b0     26.9108 28.9716 30.0698 31.166 33.2015
## b1      0.4218  0.4863  0.5211  0.556  0.6222
## sigma2 20.3575 24.1980 26.6129 29.370 35.6623

autocorr.diag(samples)

##                b0         b1      sigma2
## Lag 0  1.00000000 1.00000000 1.000000000
## Lag 1  0.89674469 0.89729301 0.018382972
## Lag 5  0.58881031 0.58970557 0.002997731
## Lag 10 0.34899078 0.35034313 0.004552666
## Lag 50 0.03019034 0.03073331 0.004531337

autocorr.plot(samples[[1]])

การอนุมานเชิงสถิติแบบเบส์

จากผลการวินิจฉัย MCMC ในข้างต้น เราทำการประมวลผลใหม่ดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

#dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
#initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model

update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 30000,
                      thin=15)
MCMCsummary(samples)

##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     29.9358012 1.61122865 26.7793446 29.9515443 33.0316408    1  4086
## b1      0.5251646 0.05136625  0.4260559  0.5245155  0.6249333    1  4067
## sigma2 26.9815883 3.93464366 20.3197188 26.6338984 35.8743312    1  6000

MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

autocorr.diag(samples)

##                   b0           b1       sigma2
## Lag 0    1.000000000 1.0000000000  1.000000000
## Lag 15   0.189732427 0.1917814786 -0.003137443
## Lag 75   0.026335394 0.0191889064  0.007447446
## Lag 150  0.006470748 0.0103972916 -0.007191825
## Lag 750 -0.003792101 0.0003762561  0.025105044

autocorr.plot(samples[[1]])

summary(samples)

## 
## Iterations = 2015:32000
## Thinning interval = 15 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## b0     29.9358 1.61123 0.0208009      0.0252067
## b1      0.5252 0.05137 0.0006631      0.0008054
## sigma2 26.9816 3.93464 0.0507960      0.0507991
## 
## 2. Quantiles for each variable:
## 
##           2.5%   25%     50%    75%   97.5%
## b0     26.7793 28.84 29.9515 31.033 33.0316
## b1      0.4261  0.49  0.5245  0.561  0.6249
## sigma2 20.3197 24.21 26.6339 29.332 35.8743

#library(googledrive)
setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel")
#drive_download("https://drive.google.com/file/d/0B3lh4V2Mrl14SF94US1YTkNVaHM/view?resourcekey=0-bYI7GCSuypi3a5aFt0O-qg")
source("DBDA2E-utilities.R") # (Kruschke, 2015)

## 
## *********************************************************************
## Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
## A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
## *********************************************************************

plotPost(samples[,"b1"], compVal=0, cenTend="median", ROPE=c(0,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.       0          1       0      0.1       0       0       1

plotPost(samples[,"b1"], compVal=0.6, cenTend="median", ROPE=c(0.55,0.65), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh   pLtROPE   pInROPE     pGtROPE
## Param. Val.     0.6 0.07133333    0.55     0.65 0.6831667 0.3091667 0.007666667

ตัวอย่างโยนเหรียญ (รอบที่ 2)

• จากตัวอย่างที่ 1 จะได้ว่าโมเดลของค่าสังเกตเป็นฟังก์ชันความน่าจะเป็นแบบ Bernuolli ที่มีพารามิเตอร์เป็น \(\theta\) เขียนแทนด้วย

• อย่างไรก็ตามในกรณีนี้ผู้วิเคราะห์ต้องการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบต่อเนื่องให้กับพารามิเตอร์ความลำเอียง \(\theta\) เนื่องจากพารามิเตอร์ดังกล่าวต้องมีค่าอยู่บนช่วง \([0,1]\) การแจกแจงความน่าจะเป็นที่เหมาะสมจึงเป็นการแจกแจงแบบบีต้า (beta distribution) ที่มีพารามิเตอร์ตำแหน่ง \(a\) และพารามิเตอร์สเกล \(b\) ดังนั้น

พารามิเตอร์ทั้งสองตัวของการแจกแจงแบบบีต้าทำให้รูปทรงการแจกแจงมีความแตกต่างกัน ซึ่งสามารถใช้สะท้อนความเชื่อเกี่ยวกับความลำเอียงของเหรียญที่แตกต่างกันได้เป็นอย่างดี รูปด้านล่างแสดงตัวอย่างโค้งความหนาแน่นของการแจกแจงแบบบีต้า เมื่อกำหนดพารามิเตอร์ต่าง ๆ

Kruschke (2012)

จากการกำหนดในข้างต้นจะสามารถเขียน syntax ของโมเดลในโปรแรกม JAGS ได้ดังนี้

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~beta(a,b)

# deterministic part
  a<-3
  b<-2
}

การประมวลผลโมเดลข้างต้นด้วยโปรแกรม JAGS บน R จำเป็นต้องมีการ wrap up คำสั่งข้างต้นให้เป็น text file เพื่อส่งออกไปยังโปรแกรม R โดยเขียนคำสั่งดังนี้

coin_model="

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~dbeta(a,b)

# deterministic part
  a<-3
  b<-2
}

"
writeLines(coin_model, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt")

ขั้นตอนถัดมาคือการประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) โดยใช้เทคนิค MCMC ผ่านโปรแกรม JAGS ซึ่งสามารถสั่งงานได้ผ่านฟังก์ชัน jag.model() โดยฟังก์ชันนี้มีอาร์กิวเมนท์ที่สำคัญที่เกี่ยวข้องกับการระบุคุณลักษณะของลูกโซ่มาร์คอฟที่ต้องการสร้างขึ้น ได้แก่

• ข้อมูลค่าสังเกต ในที่นี้สมมุติว่าทำการทดลอง 20 ครั้งได้ผลเป็นดังนี้

y<-rbinom(20,1,0.65)
data.list<-list(y=y, n=20)

• n.chains

library(rjags)
fit<-jags.model(file="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt", data=data.list,
                n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 20
##    Unobserved stochastic nodes: 1
##    Total graph size: 24
## 
## Initializing model

samples<-coda.samples(fit,variable.names=c("theta"),n.iter=5000)

plot(samples)

autocorr.diag(samples)

##              theta
## Lag 0   1.00000000
## Lag 1   0.01582687
## Lag 5  -0.01269683
## Lag 10  0.01048772
## Lag 50 -0.01040355

summary(samples)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.6796919      0.0908805      0.0007420      0.0007499 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.4931 0.6189 0.6847 0.7465 0.8412

เหรียญเที่ยงตรงหรือไม่?? —> กำหนด ROPE = [0.48,0.52]

plotPost(samples[,"theta"], compVal=0.5, cenTend="median", ROPE=c(0.48,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 14532.27 0.6796919 0.6847037 0.7075729    0.95 0.5029832 0.8492614
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.     0.5     0.9702    0.48      0.1  0.0186       0       1

Null model (random-effect ANOVA)

random-intercept model

random-coefficient model